#-----------------------------------------------------------------------
# Title: Statistical models, tables, figures for:
# "International Constraints, Political Turnover, and Voting Consistency 
# in the United Nations General Assembly", accepted for publication in Foreign 
# Policy Analysis
# Authors: Matthew DiLorenzo and Bryan Rooney
# Contact: mdiloren@odu.edu
# Institution: Old Dominion University
# Date: Last updated 04/20/2020
# Description: This script replicates all of the tables and figures
# in the main text and Online Appendix of the article.
# ----------------------------------------------------------------------

## Set working directory. To replicate, set the working directory to the
## location of "replication-package" on your machine.
setwd("~/Dropbox/article-manuscripts/conditioning-effects-fpa/final-files/replication-package/")


## Uncomment if package installation is needed
# install.packages("tidyverse")
# install.packages("stargazer")
# install.packages("xtable")
# install.packages("MASS")
# install.packages("gridExtra")

library(tidyverse)
library(stargazer)
library(xtable)

## Load custom function interaction confidence intervals and factor conversion
source("aux/extra-functions.R")

## Load data
dat <- read.csv("data/conditioning-effects-revised.csv", 
                stringsAsFactors = FALSE)

## Make GDP per capita variable
dat$gdp_per_capita <- log((dat$realgdp / dat$pop) + 1)



## Interaction models + figures --------------------------------------------

# Model 1: Number of defensive allies
sfit_1 <- lm(logvotech ~
               solsch * n_defenders +
               othldrtrans * n_defenders + 
               demboth +
               gdp_per_capita +
               ln_population +  
               regtrans + 
               CWend +
               USally +
               USSRally +
               factor(ccode) - 1, 
             data = dat[dat$interim == 0, ])

## Get upper bound of number of defensive allies to consider
n_def_max <- round(
  mean(sfit_1$model$n_defenders) +
    2 * sd(sfit_1$model$n_defenders), 0)

## Get unique values of defensive ally count
def_values <- unique((sfit_1$model$n_defenders))

## Extract model data frame
plot_data <- sfit_1$model

## Placeholder for coefficients / SEs
sfit_1_plot <- c()

for(i in 1:length(def_values)){
  sfit_1_plot <- rbind(sfit_1_plot,
                       mcce_ci_calc("solsch", 
                                    "n_defenders", 
                                    def_values[i], 
                                    sfit_1))
}

## Convert to data.frame, add defensive ally values
sfit_1_plot <- as.data.frame(sfit_1_plot)
sfit_1_plot$n_defenders <- def_values

## Merge with plot data
plot_data <- left_join(plot_data, sfit_1_plot)

## Make plot
allies_figure <- ggplot(plot_data,
                        aes(x = n_defenders, 
                            y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, n_defenders < n_def_max + 1),
              aes(x = n_defenders, 
                  y = beta,
                  ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) + 
  labs(x = "Number of defensive allies",  
       y = "Marginal effect of SOLS change",
       title = "Defensive alliances",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        "defensive allies to ",
                        "2 standard deviations above mean."))


## Save plot (uncomment to save individual plot)
# ggsave(plot = allies_figure,
#        "figures/sols-defense.pdf",
#        height = 4,
#        width = 8)

#------------------------------------------------------------
# MODEL 2: Number of rivals
#------------------------------------------------------------
sfit_2 <- lm(logvotech ~
               solsch * n_rivals +
               othldrtrans * n_rivals + 
               demboth +
               gdp_per_capita +
               ln_population +
               regtrans + 
               CWend +
               USally +
               USSRally +
               factor(ccode) - 1, 
             data = dat[dat$interim == 0, ])

n_rivals_max <- round(mean(sfit_2$model$n_rivals) +
                        2 * sd(sfit_2$model$n_rivals), 0)

## Get unique values of rivals count
riv_values <- unique((sfit_2$model$n_rivals))

## Extract model data frame
plot_data <- sfit_2$model

## Placeholder for coefficients / SEs
sfit_2_plot <- c()

for(i in 1:length(riv_values)){
  sfit_2_plot <- rbind(sfit_2_plot,
                       mcce_ci_calc("solsch", 
                                    "n_rivals", 
                                    riv_values[i], 
                                    sfit_2))
}

## Convert to data.frame, add defensive ally values
sfit_2_plot <- as.data.frame(sfit_2_plot)
sfit_2_plot$n_rivals <- riv_values

## Merge with plot data
plot_data <- left_join(plot_data, sfit_2_plot)

## Make plot
rivals_figure <- ggplot(plot_data,
                        aes(x = n_rivals, 
                            y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, n_rivals < n_rivals_max + 1),
              aes(x = n_rivals, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Number of rivals", 
       y = "Marginal effect of SOLS change",
       title = "Rivalries",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        "number of rivals to ",
                        "2 standard deviations above mean (3)."))

## Save plot (uncomment to save individual plot)
# ggsave(plot = rivals_figure,
#        "figures/sols-rivals.pdf",
#        height = 4,
#        width = 8)

#-----------------------------------------------------------------------
# MODEL 3: State capacity
#-----------------------------------------------------------------------
capacity_fit_1 <- lm(logvotech ~
                       solsch * state_capacity +
                       othldrtrans * state_capacity+ 
                       demboth +
                       gdp_per_capita +
                       ln_population +
                       regtrans + 
                       CWend +
                       USally +
                       USSRally +
                       factor(ccode) - 1, 
                     data = dat[dat$interim == 0, ])

capacity_max <- round(mean(capacity_fit_1$model$state_capacity) +
                        2*sd(capacity_fit_1$model$state_capacity), 0)

capacity_min <- round(mean(capacity_fit_1$model$state_capacity) -
                        2*sd(capacity_fit_1$model$state_capacity), 0)

## Get unique values of capacity
capacity_values <- unique((capacity_fit_1$model$state_capacity))

## Extract model data frame
plot_data <- capacity_fit_1$model

## Placeholder for coefficients / SEs
capacity_fit_1_plot <- c()

## Get marginal effects (takes a couple minutes to run)
for(i in 1:length(capacity_values)){
  capacity_fit_1_plot <- rbind(capacity_fit_1_plot,
                               mcce_ci_calc("solsch", 
                                            "state_capacity", 
                                            capacity_values[i], 
                                            capacity_fit_1))
}

## Convert to data.frame, add defensive ally values
capacity_fit_1_plot <- as.data.frame(capacity_fit_1_plot)
capacity_fit_1_plot$state_capacity <- capacity_values


## Merge with plot data
plot_data <- left_join(plot_data, capacity_fit_1_plot)

## Make plot
capacity_figure <- ggplot(plot_data, 
                          aes(x = state_capacity, 
                              y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, 
                            (state_capacity <= capacity_max) &
                              (state_capacity >= 0)),
              aes(x = state_capacity, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "State capacity", 
       y = "Marginal effect of SOLS change",
       title = "State capacity",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfor ",
                        "2 standard deviations above/below mean."))

## Save plot (uncomment to save individual plot)
# ggsave(plot = capacity_figure,
#        "figures/sols-capacity.pdf",
#        height = 4,
#        width = 8)



## Save all as one grid for final version of Figure 1

library(gridExtra)

ggsave(
  plot = grid.arrange(allies_figure, rivals_figure, capacity_figure, ncol = 2),
  filename = "figures/fig-1.png",
  width = 10, height = 8
)




#------------------------------------------------------------
# MODEL 4: Total IGO memberships
#------------------------------------------------------------
## IGO model
igo_fit_1 <- lm(logvotech ~
                  solsch * total_igos_impute +
                  othldrtrans * total_igos_impute + 
                  demboth +
                  gdp_per_capita +
                  ln_population +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])

## Get upper bound for + 2 SDs
n_igos_max <- round(mean(igo_fit_1$model$total_igos_impute) +
                      2 * sd(igo_fit_1$model$total_igos_impute), 0)

## Get unique values of IGOs count
igo_values <- unique((igo_fit_1$model$total_igos_impute))

## Extract model data frame
plot_data <- igo_fit_1$model


## Placeholder for coefficients / SEs
igo_fit_1_plot <- c()

for(i in 1:length(igo_values)){
  igo_fit_1_plot <- rbind(igo_fit_1_plot,
                          mcce_ci_calc("solsch", 
                                       "total_igos_impute", 
                                       igo_values[i], 
                                       igo_fit_1))
}

## Convert to data.frame, add defensive ally values
igo_fit_1_plot <- as.data.frame(igo_fit_1_plot)
igo_fit_1_plot$total_igos_impute <- igo_values

## Merge with plot data
plot_data <- left_join(plot_data, igo_fit_1_plot)

## Make plot
igos_figure <- ggplot(plot_data,
                      aes(x = total_igos_impute, 
                          y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, total_igos_impute < n_igos_max + 1),
              aes(x = total_igos_impute, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Number of IGO memberships", 
       y = "Marginal effect of SOLS change",
       title = "IGO Memberships",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        "number of IGOs to ",
                        "2 standard deviations above mean (89)."))

## Save plot (uncomment to save individual plot)
# ggsave(plot = igos_figure,
#        filename = "figures/sols-igos.pdf",
#        height = 4,
#        width = 8)



#-----------------------------------------------------------------------
# MODEL 5: Trade dependence
#-----------------------------------------------------------------------
trade_fit_1 <- lm(logvotech ~
                    solsch * trade_pct_gdp +
                    othldrtrans * trade_pct_gdp+ 
                    demboth +
                    gdp_per_capita +
                    ln_population +
                    regtrans + 
                    CWend +
                    USally +
                    USSRally +
                    factor(ccode) - 1, 
                  data = dat[dat$interim == 0, ])

trade_max <- round(mean(trade_fit_1$model$trade_pct_gdp) +
                     2*sd(trade_fit_1$model$trade_pct_gdp), 0)

trade_min <- round(mean(trade_fit_1$model$trade_pct_gdp) -
                     2*sd(trade_fit_1$model$trade_pct_gdp), 0)

## Get unique values of capacity
trade_values <- unique((trade_fit_1$model$trade_pct_gdp))

## Extract model data frame
plot_data <- trade_fit_1$model

## Placeholder for coefficients / SEs
trade_fit_1_plot <- c()

## Get marginal effects (takes a couple minutes to run)
for(i in 1:length(trade_values)){
  trade_fit_1_plot <- rbind(trade_fit_1_plot,
                            mcce_ci_calc("solsch", 
                                         "trade_pct_gdp", 
                                         trade_values[i], 
                                         trade_fit_1))
}

## Convert to data.frame, add defensive ally values
trade_fit_1_plot <- as.data.frame(trade_fit_1_plot)
trade_fit_1_plot$trade_pct_gdp <- trade_values

## Merge with plot data
plot_data <- left_join(plot_data, trade_fit_1_plot)

## Make plot
trade_figure <- ggplot(plot_data, 
                       aes(x = trade_pct_gdp, 
                           y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, 
                            (trade_pct_gdp <= trade_max) &
                              (trade_pct_gdp >= 0)),
              aes(x = trade_pct_gdp, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Trade as a percent of GDP", 
       y = "Marginal effect of SOLS change",
       title = "Trade dependence",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfor ",
                        "2 standard deviations above/below mean."))

## Save plot (uncomment to save individual plot)
# ggsave(plot = trade_figure,
#        filename = "figures/sols-trade.pdf",
#        height = 4,
#        width = 8)


## Save all as one grid for final version of Figure 2
ggsave(
  plot = grid.arrange(igos_figure, trade_figure, ncol = 2),
  filename = "figures/fig-2.png",
  width = 10, height = 4
)

#-----------------------------------------------------------------------
# MODEL 6: Aid
#-----------------------------------------------------------------------
aid_fit_1 <- lm(logvotech ~
                  solsch * ln_total_budget_support +
                  othldrtrans * ln_total_budget_support + 
                  demboth +
                  gdp_per_capita +
                  ln_population +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])


## Get upper bound for + 2 SDs
aid_max <- round(mean(aid_fit_1$model$ln_total_budget_support) +
                   2 * sd(aid_fit_1$model$ln_total_budget_support), 0)

aid_min <- round(mean(aid_fit_1$model$ln_total_budget_support) -
                   2 * sd(aid_fit_1$model$ln_total_budget_support), 0)

## Minus 2 SDs is -11, replacing as 0
aid_min <- 0

## Get unique values of aid
aid_values <- unique((aid_fit_1$model$ln_total_budget_support))

## Extract model data frame
plot_data <- aid_fit_1$model


## Placeholder for coefficients / SEs
aid_fit_1_plot <- c()

## Get marginal effects (takes a minute to run)
for(i in 1:length(aid_values)){
  aid_fit_1_plot <- rbind(aid_fit_1_plot,
                          mcce_ci_calc("solsch", 
                                       "ln_total_budget_support", 
                                       aid_values[i], 
                                       aid_fit_1))
}

## Convert to data.frame, add defensive ally values
aid_fit_1_plot <- as.data.frame(aid_fit_1_plot)
aid_fit_1_plot$ln_total_budget_support <- aid_values


## Merge with plot data
plot_data <- left_join(plot_data, aid_fit_1_plot)


## Make plot
aid_figure <- ggplot(plot_data,
                     aes(x = ln_total_budget_support, 
                         y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, 
                            (ln_total_budget_support < aid_max + 1) & 
                              (ln_total_budget_support > aid_min)),
              aes(x = ln_total_budget_support, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Total budget support (log)", 
       y = "Marginal effect of SOLS change",
       title = "Foreign aid",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        "aid to 2 standard deviations above mean."))

## Save plot (uncomment to save individual plot)
# ggsave(plot = aid_figure,
#        filename = "figures/sols-aid.pdf",
#        height = 4,
#        width = 8)

#-----------------------------------------------------------------------
# MODEL 7: FDI
#-----------------------------------------------------------------------
fdi_fit_1 <- lm(logvotech ~
                  solsch * fdi_pct_gdp +
                  othldrtrans * fdi_pct_gdp + 
                  demboth +
                  gdp_per_capita +
                  ln_population +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])


fdi_max <- round(mean(fdi_fit_1$model$fdi_pct_gdp) +
                   2*sd(fdi_fit_1$model$fdi_pct_gdp), 0)

fdi_min <- round(mean(fdi_fit_1$model$fdi_pct_gdp) -
                   2*sd(fdi_fit_1$model$fdi_pct_gdp), 0)

## Get unique values of fdi
fdi_values <- unique((fdi_fit_1$model$fdi_pct_gdp))

## Extract model data frame
plot_data <- fdi_fit_1$model


## Placeholder for coefficients / SEs
fdi_fit_1_plot <- c()

## Get marginal effects (takes a couple minutes to run)
for(i in 1:length(fdi_values)){
  fdi_fit_1_plot <- rbind(fdi_fit_1_plot,
                          mcce_ci_calc("solsch", 
                                       "fdi_pct_gdp", 
                                       fdi_values[i], 
                                       fdi_fit_1))
}

## Convert to data.frame, add defensive ally values
fdi_fit_1_plot <- as.data.frame(fdi_fit_1_plot)
fdi_fit_1_plot$fdi_pct_gdp <- fdi_values


## Merge with plot data
plot_data <- left_join(plot_data, fdi_fit_1_plot)


## Make plot
fdi_figure <- ggplot(subset(plot_data, 
                            (fdi_pct_gdp <= 100) & 
                              (fdi_pct_gdp > 0)), ## Using only < 100 for visual
                     aes(x = fdi_pct_gdp, 
                         y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, 
                            (fdi_pct_gdp < fdi_max + 1) &
                              (fdi_pct_gdp > 0)),
              aes(x = fdi_pct_gdp, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Foreign direct investment as a percent of GDP", 
       y = "Marginal effect of SOLS change",
       title = "FDI",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom ",
                        "2 standard deviations above/below mean."))

## Save plot (uncomment to save individual plot)
# ggsave(plot = fdi_figure,
#        filename = "figures/sols-fdi.pdf",
#        height = 4,
#        width = 8)

ggsave(
  plot = grid.arrange(aid_figure, fdi_figure, ncol = 2),
  filename = "figures/fig-3.png",
  width = 10, height = 4
)


#-----------------------------------------------------------------------
# Export models to table
#-----------------------------------------------------------------------

covariates <- c(
  "SOLS Change",
  "No. Defenders",
  "No. Rivals",
  "State Capacity",
  "Total IGO Memberships",
  "Trade (pct. GDP)",
  "Aid (log)",
  "FDI/GDP",
  "Democracy",
  "GDP per capita",
  "Population (log)",
  "Regime Transition",
  "Cold War End",
  "US Ally",
  "USSR Ally",
  "SOLS Change $\\times$ No. Defenders",
  "SOLS Change $\\times$ No. Rivals",
  "SOLS Change $\\times$ State Capacity",
  "SOLS Change $\\times$ No. IGOs",
  "SOLS Change $\\times$ Trade (pct. GDP)",
  "SOLS Change $\\times$ Aid",
  "SOLS Change $\\times$ FDI/GDP"
)


getwd()
stargazer(sfit_1, sfit_2, capacity_fit_1, 
          igo_fit_1, trade_fit_1,
          aid_fit_1, fdi_fit_1,
          dep.var.labels = c("Change in UNGA Ideal Point"),
          omit = c("ccode", "othldrtrans"),
          font.size = "tiny",
          covariate.labels = covariates,
          label = "interaction-table",
          title = c(
            paste0("International context, domestic turnover, ",
                   "and foreign policy change")),
          column.labels = c("Defense Pacts",
                            "Rivalries",
                            "Capacity",
                            "IGO Memberships",
                            "Trade",
                            "Aid",
                            "FDI / GDP"),
          single.row = FALSE,
          no.space = TRUE,
          notes.align = "l",
          notes.label = "",
          df = FALSE,
          notes = c(
            "Two-tailed tests. Estimated standard errors in parentheses.", 
            "OLS estimates. Country dummies included in all models."),
          out = c("tables/interaction-table.tex", 
                  "tables/interaction-table.html"))



#-----------------------------------------------------------------------
# Main text: Three-way interaction models + figures -- regime type
#-----------------------------------------------------------------------

# Model 1: Number of defensive allies
s3fit_1 <- lm(logvotech ~
                solsch * n_defenders * demboth +
                othldrtrans * n_defenders * demboth + 
                gdp_per_capita +
                ln_population +
                regtrans + 
                CWend +
                USally +
                USSRally +
                factor(ccode) - 1, 
              data = dat[dat$interim == 0, ])

## Get upper bound of number of defensive allies to consider
n_def_max <- round(mean(s3fit_1$model$n_defenders) +
                     2 * sd(s3fit_1$model$n_defenders), 0)

defenders_range <- c(0, n_def_max)

# Model 2: Number of rivals

s3fit_2 <- lm(logvotech ~
                solsch * n_rivals * demboth +
                othldrtrans * n_rivals * demboth + 
                gdp_per_capita +
                ln_population +
                regtrans + 
                CWend +
                USally +
                USSRally +
                factor(ccode) - 1, 
              data = dat[dat$interim == 0, ])

n_rivals_max <- round(mean(s3fit_2$model$n_rivals) +
                        2 * sd(s3fit_2$model$n_rivals), 0)

rivals_range <- c(0, n_rivals_max)

## State capacity -------------------------------------------------------------

capacity_3fit_1 <- lm(logvotech ~
                        solsch * state_capacity * demboth +
                        othldrtrans * state_capacity * demboth + 
                        gdp_per_capita +
                        ln_population +
                        regtrans + 
                        CWend +
                        USally +
                        USSRally +
                        factor(ccode) - 1, 
                      data = dat[dat$interim == 0, ])


capacity_max <- round(mean(capacity_3fit_1$model$state_capacity) +
                        2*sd(capacity_3fit_1$model$state_capacity), 0)

capacity_min <- round(mean(capacity_3fit_1$model$state_capacity) -
                        2*sd(capacity_3fit_1$model$state_capacity), 0)

capacity_range <- c(capacity_min, capacity_max)

## IGO model ----------------------------------------------------------------
igo_3fit_1 <- lm(logvotech ~
                   solsch * total_igos_impute * demboth +
                   othldrtrans * total_igos_impute * demboth + 
                   gdp_per_capita +
                   ln_population +
                   regtrans + 
                   CWend +
                   USally +
                   USSRally +
                   factor(ccode) - 1, 
                 data = dat[dat$interim == 0, ])

## Get upper bound for + 2 SDs
n_igos_max <- round(mean(igo_3fit_1$model$total_igos_impute) +
                      2 * sd(igo_3fit_1$model$total_igos_impute), 0)


igos_range <- c(0, n_igos_max)

## Political economy models -------------------------------------------------
aid_3fit_1 <- lm(logvotech ~
                   solsch * ln_total_budget_support * demboth +
                   othldrtrans * ln_total_budget_support * demboth + 
                   gdp_per_capita +
                   ln_population +
                   regtrans + 
                   CWend +
                   USally +
                   USSRally +
                   factor(ccode) - 1, 
                 data = dat[dat$interim == 0, ])


## Get upper bound for + 2 SDs
aid_max <- round(mean(aid_3fit_1$model$ln_total_budget_support) +
                   2 * sd(aid_3fit_1$model$ln_total_budget_support), 0)

aid_min <- round(mean(aid_3fit_1$model$ln_total_budget_support) -
                   2 * sd(aid_3fit_1$model$ln_total_budget_support), 0)

## Minus 2 SDs is -11, replacing as 0
aid_min <- 0

aid_range <- c(aid_min, aid_max)



## Foreign direct investment ------------------------------------------------- 
fdi_3fit_1 <- lm(logvotech ~
                   solsch * fdi_pct_gdp * demboth +
                   othldrtrans * fdi_pct_gdp * demboth + 
                   gdp_per_capita +
                   ln_population +
                   regtrans + 
                   CWend +
                   USally +
                   USSRally +
                   factor(ccode) - 1, 
                 data = dat[dat$interim == 0, ])


fdi_max <- round(mean(fdi_3fit_1$model$fdi_pct_gdp) +
                   2*sd(fdi_3fit_1$model$fdi_pct_gdp), 0)

fdi_min <- round(mean(fdi_3fit_1$model$fdi_pct_gdp) -
                   2*sd(fdi_3fit_1$model$fdi_pct_gdp), 0)

fdi_range <- c(0, fdi_max)


## Trade dependence -----------------------------------------------------------

trade_3fit_1 <- lm(logvotech ~
                     solsch * trade_pct_gdp * demboth +
                     othldrtrans * trade_pct_gdp * demboth + 
                     gdp_per_capita +
                     ln_population +
                     regtrans + 
                     CWend +
                     USally +
                     USSRally +
                     factor(ccode) - 1, 
                   data = dat[dat$interim == 0, ])


trade_max <- round(mean(trade_3fit_1$model$trade_pct_gdp) +
                     2*sd(trade_3fit_1$model$trade_pct_gdp), 0)

trade_min <- round(mean(trade_3fit_1$model$trade_pct_gdp) -
                     2*sd(trade_3fit_1$model$trade_pct_gdp), 0)

trade_range <- c(trade_min, trade_max)




#-----------------------------------------------------------------------
# Figures
#-----------------------------------------------------------------------


ranges <- list(defenders_range,
               rivals_range,
               capacity_range,
               igos_range,
               trade_range,
               aid_range,
               fdi_range)

fits <- list(s3fit_1,
             s3fit_2,
             capacity_3fit_1,
             igo_3fit_1,
             trade_3fit_1,
             aid_3fit_1,
             fdi_3fit_1)


conditioning_values <- list(
  unique(s3fit_1$model$n_defenders),
  unique(s3fit_2$model$n_rivals),
  unique(capacity_3fit_1$model$state_capacity),
  unique(igo_3fit_1$model$total_igos_impute),
  unique(trade_3fit_1$model$trade_pct_gdp),
  unique(aid_3fit_1$model$ln_total_budget_support),
  unique(fdi_3fit_1$model$fdi_pct_gdp)
)

titles <- c("Number of defensive allies",
            "Number of rivals",
            "State capacity",
            "Number of IGO memberships",
            "Trade",
            "Aid dependence",
            "FDI")

## Simulate betas with mean = coefs from model and model var-cov matrix
library(MASS)

separate_figures <- list()

for(i in 1:length(fits)){
  
  sim_beta <- mvrnorm(n = 1000, 
                      mu = coef(fits[[i]]),
                      Sigma = vcov(fits[[i]]))
  
  ## Keep necessary columns
  sim_beta <- sim_beta %>%
    as.data.frame() %>%
    dplyr::select(contains("solsch"))
  
  ## Place holder for each data set
  calc_data <- list()
  
  for(j in 1:nrow(sim_beta)){
    
    betas <- sim_beta[j, ] %>% unlist() %>% c()
    
    ## Set up data set to look at all combos of condition / democracy
    calc_data[[j]] <- cbind.data.frame(
      condition = rep(conditioning_values[[i]], 2),
      Democracy = c(rep(0, length(conditioning_values[[i]])),
                    rep(1, length(conditioning_values[[i]]))),
      b_1 = betas[1],
      b_2 = betas[2],
      b_3 = betas[3],
      b_4 = betas[4],
      row.names = NULL
    ) %>%
      mutate(
        beta_sols = b_1 + 
          (b_2 * condition) + 
          (b_3 * Democracy) + 
          (b_4 * condition * Democracy)
      ) %>%
      dplyr::select(condition, Democracy, beta_sols) %>%
      mutate(iteration = j)
    
  }
  
  ## Bind all together
  calc_data <- do.call(rbind, calc_data) %>%
    as.data.frame()
  
  ## Calculate mean and lower/upper quantile
  calc_data <- calc_data %>%
    group_by(condition, Democracy) %>%
    summarise(mean_beta = mean(beta_sols),
              beta_lower = quantile(beta_sols, 0.025),
              beta_upper = quantile(beta_sols, 0.975)) %>%
    mutate(Democracy = recode(Democracy,
                              `0` = "Non-democracy",
                              `1` = "Democracy")) %>%
    ungroup() %>%
    rename(`Regime Type` = Democracy) 
  
  ## Plot
  p <- ggplot(
    calc_data %>%
      subset(condition <= ranges[[i]][2] & condition >= ranges[[i]][1]),
    aes(x = condition, y = mean_beta,
        ymin = beta_lower,
        ymax = beta_upper,
        group = `Regime Type`,
        fill = `Regime Type`)) +
    geom_line(stat = "identity") +
    geom_hline(yintercept = 0, lty = "dashed") +
    geom_ribbon(alpha = 0.35) +
    scale_fill_grey(start = 0.8, end = 0.2) +
    theme_plot_plain() +
    labs(x = "Value of conditioning variable",
         y = "Estimated marginal effect of a SOLS change",
         title = titles[i]) 
  
  ## Uncomment to save each figure separately
  # the_file <- paste0("figures/",
  #                    titles[i] %>% tolower() %>% gsub(" ", "-", .),
  #                    "-dem-vs-nondem",
  #                    ".jpeg")
  # 
  # ## Save
  # ggsave(p,
  #        file = the_file)
  
  separate_figures[[i]] <- p
  
}


## Make full plot
library(gridExtra)

ggsave(
  plot = grid.arrange(
    separate_figures[[1]],
    separate_figures[[2]],
    separate_figures[[3]],
    separate_figures[[4]],
    separate_figures[[5]],
    separate_figures[[6]],
    separate_figures[[7]], 
    ncol = 2
  ),
  filename = "figures/fig-4.png",
  width = 12, height = 16
)


## Export model output (copied into .tex file for Online Appendix)
stargazer(
  fits,
  title = "Output for three-way interaction models",
  omit = c("ccode", "othldr"),
  no.space = TRUE,
  omit.stat = c("f", "ser"),
  font.size = "tiny",
  out = c("appendix-tables/three-way-interaction-output.tex",
          "appendix-tables/three-way-interaction-output.html")
)


#-----------------------------------------------------------------------
# Main text: Figure 5 -- three-way interaction between SOLS, conditions, GDP
#-----------------------------------------------------------------------

# Model 1: Number of defensive allies
sfit_3waygdp_1 <- lm(logvotech ~
                       solsch * n_defenders * gdp_per_capita +
                       othldrtrans * n_defenders * gdp_per_capita + 
                       demboth + 
                       ln_population + 
                       regtrans + 
                       CWend +
                       USally +
                       USSRally +
                       factor(ccode) - 1, 
                     data = dat[dat$interim == 0, ])

## Get upper bound of number of defensive allies to consider
n_def_max <- round(mean(sfit_3waygdp_1$model$n_defenders) +
                     2 * sd(sfit_3waygdp_1$model$n_defenders), 0)

defenders_range <- c(0, n_def_max)

# Model 2: Number of rivals

sfit_3waygdp_2 <- lm(logvotech ~
                       solsch * n_rivals * gdp_per_capita +
                       othldrtrans * n_rivals * gdp_per_capita + 
                       demboth + 
                       ln_population +
                       regtrans + 
                       CWend +
                       USally +
                       USSRally +
                       factor(ccode) - 1, 
                     data = dat[dat$interim == 0, ])

n_rivals_max <- round(mean(sfit_3waygdp_2$model$n_rivals) +
                        2 * sd(sfit_3waygdp_2$model$n_rivals), 0)

rivals_range <- c(0, n_rivals_max)

## State capacity -------------------------------------------------------------

capacity_fit_3waygdp_1 <- lm(logvotech ~
                               solsch * state_capacity * gdp_per_capita +
                               othldrtrans * state_capacity * gdp_per_capita + 
                               demboth + 
                               ln_population +
                               regtrans + 
                               CWend +
                               USally +
                               USSRally +
                               factor(ccode) - 1, 
                             data = dat[dat$interim == 0, ])


capacity_max <- round(mean(capacity_fit_3waygdp_1$model$state_capacity) +
                        2*sd(capacity_fit_3waygdp_1$model$state_capacity), 0)

capacity_min <- round(mean(capacity_fit_3waygdp_1$model$state_capacity) -
                        2*sd(capacity_fit_3waygdp_1$model$state_capacity), 0)

capacity_range <- c(capacity_min, capacity_max)

## IGO model ----------------------------------------------------------------
igo_fit_3waygdp_1 <- lm(logvotech ~
                          solsch * total_igos_impute * gdp_per_capita +
                          othldrtrans * total_igos_impute * gdp_per_capita + 
                          demboth + 
                          ln_population +
                          regtrans + 
                          CWend +
                          USally +
                          USSRally +
                          factor(ccode) - 1, 
                        data = dat[dat$interim == 0, ])

## Get upper bound for + 2 SDs
n_igos_max <- round(mean(igo_fit_3waygdp_1$model$total_igos_impute) +
                      2 * sd(igo_fit_3waygdp_1$model$total_igos_impute), 0)


igos_range <- c(0, n_igos_max)

## Political economy models -------------------------------------------------
aid_fit_3waygdp_1 <- lm(logvotech ~
                          solsch * ln_total_budget_support * gdp_per_capita +
                          othldrtrans * ln_total_budget_support * gdp_per_capita + 
                          demboth + 
                          ln_population + 
                          regtrans + 
                          CWend +
                          USally +
                          USSRally +
                          factor(ccode) - 1, 
                        data = dat[dat$interim == 0, ])


## Get upper bound for + 2 SDs
aid_max <- round(mean(aid_fit_3waygdp_1$model$ln_total_budget_support) +
                   2 * sd(aid_fit_3waygdp_1$model$ln_total_budget_support), 0)

aid_min <- round(mean(aid_fit_3waygdp_1$model$ln_total_budget_support) -
                   2 * sd(aid_fit_3waygdp_1$model$ln_total_budget_support), 0)

## Minus 2 SDs is -11, replacing as 0
aid_min <- 0
aid_range <- c(aid_min, aid_max)



## Foreign direct investment ------------------------------------------------- 
fdi_fit_3waygdp_1 <- lm(logvotech ~
                          solsch * fdi_pct_gdp * gdp_per_capita +
                          othldrtrans * fdi_pct_gdp * gdp_per_capita + 
                          demboth + 
                          ln_population + 
                          regtrans + 
                          CWend +
                          USally +
                          USSRally +
                          factor(ccode) - 1, 
                        data = dat[dat$interim == 0, ])


fdi_max <- round(mean(fdi_fit_3waygdp_1$model$fdi_pct_gdp) +
                   2*sd(fdi_fit_3waygdp_1$model$fdi_pct_gdp), 0)

fdi_min <- round(mean(fdi_fit_3waygdp_1$model$fdi_pct_gdp) -
                   2*sd(fdi_fit_3waygdp_1$model$fdi_pct_gdp), 0)

fdi_range <- c(0, fdi_max)


## Trade dependence -----------------------------------------------------------

trade_fit_3waygdp_1 <- lm(logvotech ~
                            solsch * trade_pct_gdp * gdp_per_capita +
                            othldrtrans * trade_pct_gdp * gdp_per_capita + 
                            demboth + 
                            ln_population + 
                            regtrans + 
                            CWend +
                            USally +
                            USSRally +
                            factor(ccode) - 1, 
                          data = dat[dat$interim == 0, ])


trade_max <- round(mean(trade_fit_3waygdp_1$model$trade_pct_gdp) +
                     2*sd(trade_fit_3waygdp_1$model$trade_pct_gdp), 0)

trade_min <- round(mean(trade_fit_3waygdp_1$model$trade_pct_gdp) -
                     2*sd(trade_fit_3waygdp_1$model$trade_pct_gdp), 0)

trade_range <- c(trade_min, trade_max)




#-----------------------------------------------------------------------
# Figures
#-----------------------------------------------------------------------

ranges <- list(defenders_range,
               rivals_range,
               capacity_range,
               igos_range,
               trade_range,
               aid_range,
               fdi_range)

fits <- list(sfit_3waygdp_1,
             sfit_3waygdp_2,
             capacity_fit_3waygdp_1,
             igo_fit_3waygdp_1,
             trade_fit_3waygdp_1,
             aid_fit_3waygdp_1,
             fdi_fit_3waygdp_1)


conditioning_values <- list(
  seq(defenders_range[1], defenders_range[2], 
      length.out = max(c(20, defenders_range[2]))),
  seq(rivals_range[1], rivals_range[2], 
      length.out = max(c(20, rivals_range[2]))),
  seq(capacity_range[1], capacity_range[2], 
      length.out = max(c(20, capacity_range[2]))),
  seq(igos_range[1], igos_range[2], 
      length.out = max(c(20, igos_range[2]))),
  seq(trade_range[1], trade_range[2], 
      length.out = max(c(20, trade_range[2]))),
  seq(aid_range[1], aid_range[2], 
      length.out = max(c(20, aid_range[2]))),
  seq(fdi_range[1], fdi_range[2], 
      length.out = max(c(20, fdi_range[2])))
)


titles <- c("Number of defensive allies",
            "Number of rivals",
            "State capacity",
            "Number of IGO memberships",
            "Trade",
            "Aid dependence",
            "FDI")

gdp_upper <- round(mean(dat$gdp_per_capita, na.rm = T) +
                     2*sd(dat$gdp_per_capita, na.rm = T), 2)

gdp_range <- c(min(dat$gdp, na.rm = T), gdp_upper)

gdp_range_values <- seq(gdp_range[1], gdp_range[2], length.out = 30)


## Simulate betas with mean = coefs from model and model var-cov matrix
library(MASS)

separate_figures <- list()

for(i in 1:length(fits)){
  
  sim_beta <- mvrnorm(n = 1000, 
                      mu = coef(fits[[i]]),
                      Sigma = vcov(fits[[i]]))
  
  ## Keep necessary columns
  sim_beta <- sim_beta %>%
    as.data.frame() %>%
    dplyr::select(contains("solsch"))
  
  ## Place holder for each data set
  calc_data <- list()
  
  for(j in 1:nrow(sim_beta)){
    
    betas <- sim_beta[j, ] %>% unlist() %>% c()
    
    ## Set up data set to look at all combos of condition / democracy
    calc_data[[j]] <- cbind.data.frame(
      condition = rep(conditioning_values[[i]], length(gdp_range_values)),
      GDP = c(rep(gdp_range_values,  each = length(conditioning_values[[i]]))),
      b_1 = betas[1],
      b_2 = betas[2],
      b_3 = betas[3],
      b_4 = betas[4],
      row.names = NULL
    ) %>%
      mutate(
        beta_sols = b_1 + 
          (b_2 * condition) + 
          (b_3 * GDP) + 
          (b_4 * condition * GDP)
      ) %>%
      dplyr::select(condition, GDP, beta_sols) %>%
      mutate(iteration = j)
    
  }
  
  ## Bind all together
  calc_data <- do.call(rbind, calc_data) %>%
    as.data.frame()
  
  ## Calculate mean and lower/upper quantile
  calc_data <- calc_data %>%
    group_by(condition, GDP) %>%
    summarise(mean_beta = mean(beta_sols),
              beta_lower = quantile(beta_sols, 0.025),
              beta_upper = quantile(beta_sols, 0.975)) %>%
    ungroup() %>%
    mutate(significant = ifelse(beta_lower * beta_upper > 0, TRUE, FALSE)) %>%
    mutate(positive = ifelse(mean_beta > 0, TRUE, FALSE)) %>%
    mutate(Result = NA) %>%
    mutate(Result = replace(Result,
                            positive & significant,
                            "+, significant")) %>%
    mutate(Result = replace(Result,
                            positive & (!significant),
                            "+")) %>%
    mutate(Result = replace(Result,
                            (!positive) & (!significant),
                            "-")) %>%
    mutate(Result = replace(Result,
                            (!positive) & significant,
                            "-, significant")) %>%
    mutate(Result = factor(Result,
                           levels = rev(c("-, significant",
                                          "-",
                                          "+",
                                          "+, significant")))) 
  
  
  
  p <- ggplot(calc_data,
              aes(x = condition, y = GDP,
                  fill = Result)) +
    scale_fill_manual(values = c("#383838",
                                 "#808080",
                                 "#d3d3d3",
                                 "#f5f5f5")) +
    geom_tile() +
    theme_plot_plain() +
    labs(
      x = titles[i],
      y = "GDP per capita (log)"
    ) 
  
  
  ## Uncomment to save each figure separately
  # getwd()
  # the_file <- paste0("figures/",
  #                    titles[i] %>% tolower() %>% gsub(" ", "-", .),
  #                    "-by-gdppc",
  #                    ".jpeg")
  
  
  ## Save
  # ggsave(p,
  #        file = the_file)
  
  separate_figures[[i]] <- p
  
}


## Make full plot
library(gridExtra)

ggsave(
  plot = grid.arrange(
    separate_figures[[1]],
    separate_figures[[2]],
    separate_figures[[3]],
    separate_figures[[4]],
    separate_figures[[5]],
    separate_figures[[6]],
    separate_figures[[7]], 
    ncol = 2
  ),
  filename = "figures/fig-5.png",
  width = 12, height = 16
)


#------------------------------------------------------------
# Online Appendix: Summary stats table
#------------------------------------------------------------
summary_data <- dat %>%
  dplyr::select(logvotech,
                solsch,
                n_defenders,
                n_rivals,
                total_igos_impute,
                ln_total_budget_support,
                fdi_pct_gdp,
                trade_pct_gdp,
                state_capacity,
                othldrtrans,
                demboth,
                gdp_per_capita,
                ln_population,
                regtrans, 
                CWend,
                USally,
                USSRally) %>%
  mutate_all(funs(factor2numeric(.)))

getwd()
stargazer(summary_data,
          title = "Summary statistics for key variables",
          label = "summary-table",
          summary.stat = c("n", "mean", "sd", "min", "max"),
          out = c("appendix-tables/summary-stats.tex",
                  "appendix-tables/summary-stats.html"))

#------------------------------------------------------------
# Online Appendix: Correlation table for conditioning variables
#------------------------------------------------------------

cor_matrix <- cor(summary_data %>%
                    dplyr::select(
                      n_defenders,
                      n_rivals,
                      state_capacity,
                      total_igos_impute,
                      ln_total_budget_support,
                      fdi_pct_gdp,
                      trade_pct_gdp
                    ),
                  use = "pairwise.complete.obs")


colnames(cor_matrix) <- rownames(cor_matrix) <- c("No. Defenders",
                                                  "No. Rivals",
                                                  "State Capacity",
                                                  "Total IGO Memberships",
                                                  "Aid (log)",
                                                  "FDI/GDP",
                                                  "Trade (pct. GDP)")

stargazer(cor_matrix,
          title = "Correlations between key conditioning variables",
          label = "cor-table",
          font.size = "tiny",
          out = c("appendix-tables/correlation-matrix-for-moderators.tex",
                  "appendix-tables/correlation-matrix-for-moderators.html"))



#-----------------------------------------------------------------------
# Online appendix: Subset of non-democratic country-years
#-----------------------------------------------------------------------

full_models <- list(sfit_1, 
                    sfit_2, 
                    capacity_fit_1, 
                    igo_fit_1, 
                    trade_fit_1,
                    aid_fit_1, 
                    fdi_fit_1)

nondem_models <- list()
dem_models <- list()

for(i in 1:length(full_models)){
  
  nondem_models[[i]] <- update(
    full_models[[i]], 
    . ~ . - demboth - regtrans,
    data = dat[dat$interim == 0 & dat$dem == 0, ]
  )
  
  dem_models[[i]] <- update(
    full_models[[i]], 
    . ~ . - demboth - regtrans,
    data = dat[dat$interim == 0 & dat$dem == 1, ]
  )
  
}


nondem_covariates <- covariates[!covariates %in% c("Democracy", 
                                                   "Regime Transition")]

getwd()
stargazer(nondem_models,
          omit = c("ccode", "othldr"),
          font.size = "tiny",
          covariate.labels = nondem_covariates,
          dep.var.labels = c("Change in UNGA Ideal Point"),
          label = "interaction-table-nondem",
          title = c(
            paste0("International context, domestic turnover, ",
                   "and foreign policy change in non-democracies")),
          #dep.var.caption = "",
          column.labels = c("Defense Pacts",
                            "Rivalries",
                            "Capacity",
                            "IGO Memberships",
                            "Trade",
                            "Aid",
                            "FDI / GDP"),
          single.row = FALSE,
          no.space = TRUE,
          notes.align = "l",
          notes.label = "",
          df = FALSE,
          notes = c(
            "Two-tailed tests. Estimated standard errors in parentheses.", 
            "OLS estimates. Country dummies included in all models."),
          out = c("appendix-tables/interaction-table-nondem.tex",
                  "appendix-tables/interaction-table-nondem.html"))


getwd()
stargazer(dem_models,
          omit = c("ccode", "othldr"),
          font.size = "tiny",
          covariate.labels = nondem_covariates,
          dep.var.labels = c("Change in UNGA Ideal Point"),
          label = "interaction-table-nondem",
          title = c(
            paste0("International context, domestic turnover, ",
                   "and foreign policy change in democracies")),
          #dep.var.caption = "",
          column.labels = c("Defense Pacts",
                            "Rivalries",
                            "Capacity",
                            "IGO Memberships",
                            "Trade",
                            "Aid",
                            "FDI / GDP"),
          single.row = FALSE,
          no.space = TRUE,
          notes.align = "l",
          notes.label = "",
          df = FALSE,
          notes = c(
            "Two-tailed tests. Estimated standard errors in parentheses.", 
            "OLS estimates. Country dummies included in all models."),
          out = c("appendix-tables/interaction-table-dem.tex",
                  "appendix-tables/interaction-table-dem.html"))



#-----------------------------------------------------------------------
# Online Appendix: Controlling/interacting with GDP per capita
#-----------------------------------------------------------------------

gdp_fit_1 <- lm(logvotech ~
                  solsch * gdp_per_capita +
                  othldrtrans * gdp_per_capita + 
                  demboth +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])


gdp_per_capita_max <- round(mean(gdp_fit_1$model$gdp_per_capita) +
                              2 * sd(gdp_fit_1$model$gdp_per_capita), 0)



## Get unique values of GDP
gdp_values <- unique((gdp_fit_1$model$gdp_per_capita))

## Extract model data frame
plot_data <- gdp_fit_1$model

## Placeholder for coefficients / SEs
gdp_fit_1_plot <- c()

for(i in 1:length(gdp_values)){
  gdp_fit_1_plot <- rbind(gdp_fit_1_plot,
                          mcce_ci_calc("solsch", 
                                       "gdp_per_capita", 
                                       gdp_values[i], 
                                       gdp_fit_1))
}

## Convert to data.frame, add defensive ally values
gdp_fit_1_plot <- as.data.frame(gdp_fit_1_plot)
gdp_fit_1_plot$gdp_per_capita <- gdp_values


## Merge with plot data
plot_data <- left_join(plot_data, gdp_fit_1_plot)

## Make plot
ggplot(plot_data,
       aes(x = gdp_per_capita, 
           y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, gdp_per_capita < gdp_per_capita_max),
              aes(x = gdp_per_capita, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "GDP per capita (log)", 
       y = "Marginal effect of SOLS change",
       title = "GDP per capita (log)",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        "GDP per capita to ",
                        "2 standard deviations above mean."))

## Save plot
ggsave("appendix-figures/appx-sols-gdp-per-capita.pdf",
       height = 4,
       width = 8)

#------------------------------------------------------------
# Online Appendix: Interacting with population
#------------------------------------------------------------

pop_fit_1 <- lm(logvotech ~
                  solsch * ln_population +
                  othldrtrans * ln_population + 
                  demboth +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])

ln_population_max <- round(mean(pop_fit_1$model$ln_population) +
                             2 * sd(pop_fit_1$model$ln_population), 0)


## Get unique values of population
pop_values <- unique((pop_fit_1$model$ln_population))

## Extract model data frame
plot_data <- pop_fit_1$model


## Placeholder for coefficients / SEs
pop_fit_1_plot <- c()

for(i in 1:length(pop_values)){
  pop_fit_1_plot <- rbind(pop_fit_1_plot,
                          mcce_ci_calc("solsch", 
                                       "ln_population", 
                                       pop_values[i], 
                                       pop_fit_1))
}

## Convert to data.frame, add defensive ally values
pop_fit_1_plot <- as.data.frame(pop_fit_1_plot)
pop_fit_1_plot$ln_population <- pop_values


## Merge with plot data
plot_data <- left_join(plot_data, pop_fit_1_plot)

## Make plot
ggplot(plot_data,
       aes(x = ln_population, 
           y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, ln_population < ln_population_max),
              aes(x = ln_population, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Population (log)", 
       y = "Marginal effect of SOLS change",
       title = "Population (log)",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        "population (log) to ",
                        "2 standard deviations above mean."))

## Save plot
ggsave("appendix-figures/appx-sols-population.pdf",
       height = 4,
       width = 8)



#-----------------------------------------------------------------------
# Online Appendix: Controlling for aid as percent of GDP
#-----------------------------------------------------------------------

dat$aid_pct_gdp <- log((dat$total_budget_support / dat$realgdp) + 1)

aid_pct_gdp_fit <- lm(logvotech ~
                        solsch * aid_pct_gdp +
                        othldrtrans * aid_pct_gdp + 
                        demboth +
                        regtrans + 
                        CWend +
                        USally +
                        USSRally +
                        factor(ccode) - 1, 
                      data = dat[dat$interim == 0, ])

aid_pct_gdp_max <- round(mean(aid_pct_gdp_fit$model$aid_pct_gdp) +
                           2 * sd(aid_pct_gdp_fit$model$aid_pct_gdp), 0)

## Get unique values of aid count
aid_gdp_values <- unique((aid_pct_gdp_fit$model$aid_pct_gdp))

## Extract model data frame
plot_data <- aid_pct_gdp_fit$model


## Placeholder for coefficients / SEs
aid_pct_gdp_fit_plot <- c()

for(i in 1:length(aid_gdp_values)){
  aid_pct_gdp_fit_plot <- rbind(aid_pct_gdp_fit_plot,
                                mcce_ci_calc("solsch", 
                                             "aid_pct_gdp", 
                                             aid_gdp_values[i], 
                                             aid_pct_gdp_fit))
}

## Convert to data.frame, add defensive ally values
aid_pct_gdp_fit_plot <- as.data.frame(aid_pct_gdp_fit_plot)
aid_pct_gdp_fit_plot$aid_pct_gdp <- aid_gdp_values


## Merge with plot data
plot_data <- left_join(plot_data, aid_pct_gdp_fit_plot)


## Make plot
ggplot(plot_data,
       aes(x = aid_pct_gdp, 
           y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, aid_pct_gdp < aid_pct_gdp_max),
              aes(x = aid_pct_gdp, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Aid / GDP (log)", 
       y = "Marginal effect of SOLS change",
       title = "Aid / GDP (log)",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        "aid as percent of GDP to ",
                        "2 standard deviations above mean."))

## Save plot
ggsave("appendix-figures/appx-sols-aid-pct-gdp.pdf",
       height = 4,
       width = 8)


#-----------------------------------------------------------------------
# Online Appendix: Export real GDP and population models
#-----------------------------------------------------------------------

cap_covariates <- c(
  "SOLS Change",
  "GDP per capita (log)",
  "Population (log)",
  "Aid / GDP (log)",
  "Democracy",
  "Regime Transition",
  "Cold War End",
  "US Ally",
  "USSR Ally",
  "SOLS Change $\\times$ GDP per capita (log)",
  "SOLS Change $\\times$ Population (log)",
  "SOLS Change $\\times$ Aid / GDP (log)"
)


## Without interactions
gdp_fit_0 <- lm(logvotech ~
                  solsch +
                  gdp_per_capita +
                  othldrtrans +
                  demboth +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])

pop_fit_0 <- lm(logvotech ~
                  solsch +
                  ln_population +
                  othldrtrans +
                  demboth +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])

aid_fit_0 <- lm(logvotech ~
                  solsch +
                  aid_pct_gdp +
                  othldrtrans +
                  demboth +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])

gdp_pop_fit_0 <- lm(logvotech ~
                      solsch +
                      gdp_per_capita +
                      ln_population +
                      aid_pct_gdp +
                      othldrtrans +
                      demboth +
                      regtrans + 
                      CWend +
                      USally +
                      USSRally +
                      factor(ccode) - 1, 
                    data = dat[dat$interim == 0, ])


getwd()
stargazer(gdp_fit_0, pop_fit_0, aid_fit_0, 
          gdp_pop_fit_0, 
          gdp_fit_1, pop_fit_1, aid_pct_gdp_fit,
          dep.var.labels = c("Change in UNGA Ideal Point"),
          omit = c("ccode", "othldr"),
          font.size = "tiny",
          covariate.labels = cap_covariates,
          label = "appx-capacity-table",
          title = c(
            paste0("Alternative measures of capacity")),
          single.row = FALSE,
          no.space = TRUE,
          notes.align = "l",
          notes.label = "",
          df = FALSE,
          notes = c(
            "Two-tailed tests. Estimated standard errors in parentheses.", 
            "OLS estimates. Country dummies included in all models."),
          out = c("appendix-tables/appx-capacity-table.tex", 
                  "appendix-tables/appx-capacity-table.html"))



#-----------------------------------------------------------------------
# Online Appendix: PCA analysis
#-----------------------------------------------------------------------

pca_data <- dat %>%
  dplyr::select(
    ccode, year, 
    n_defenders,
    n_rivals,
    total_igos_impute,
    ln_total_budget_support,
    fdi_pct_gdp,
    trade_pct_gdp,
    state_capacity
  ) %>% 
  na.exclude()


pr.out <- prcomp(
  pca_data %>% 
    dplyr::select(-ccode, -year), 
  scale = TRUE)

pca_data$constraint_score <- as.numeric(
  as.matrix(pca_data %>% 
              dplyr::select(-ccode, -year)) %*% 
    as.matrix(pr.out$rotation[,1]^2)
)


pr.out$rotation
sum(pr.out$rotation[ , 2]^2)
dim(pr.out$x)

pr.var <- pr.out$sdev^2
pve <- pr.var/sum(pr.var)

jpeg(file = "appendix-figures/pca-variance-explained.jpeg")
plot(cumsum(pve),
     xlab = "Principal Component",
     ylab = "Cumulative Proportion of Variance Explained",
     ylim = c(0, 1),
     type = 'b')
dev.off()


dat <- left_join(dat, 
                 pca_data %>% 
                   dplyr::select(ccode, year, constraint_score))

pca_fit_1 <- lm(logvotech ~
                  solsch * constraint_score +
                  othldrtrans * constraint_score + 
                  demboth +
                  gdp_per_capita +
                  ln_population +
                  regtrans + 
                  CWend +
                  USally +
                  USSRally +
                  factor(ccode) - 1, 
                data = dat[dat$interim == 0, ])

constraint_score_max <- round(mean(pca_fit_1$model$constraint_score) +
                                2 * sd(pca_fit_1$model$constraint_score), 0)

## Get unique values
pca_values <- unique((pca_fit_1$model$constraint_score))

## Extract model data frame
plot_data <- pca_fit_1$model


## Placeholder for coefficients / SEs
pca_fit_1_plot <- c()

for(i in 1:length(pca_values)){
  pca_fit_1_plot <- rbind(pca_fit_1_plot,
                          mcce_ci_calc("solsch", 
                                       "constraint_score", 
                                       pca_values[i], 
                                       pca_fit_1))
}

## Convert to data.frame, add values
pca_fit_1_plot <- as.data.frame(pca_fit_1_plot)
pca_fit_1_plot$constraint_score <- pca_values


## Merge with plot data
plot_data <- left_join(plot_data, pca_fit_1_plot)


## Make plot
ggplot(plot_data,
       aes(x = constraint_score, 
           y = beta, ymin = ci.95.lower, ymax = ci.95.upper)) +
  geom_line() +
  geom_rug(sides = "b") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  geom_ribbon(data = subset(plot_data, constraint_score < constraint_score_max),
              aes(x = constraint_score, 
                  y = beta, ymin = ci.95.lower, ymax = ci.95.upper),
              alpha = .3, fill = "black") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "International constraint PCA score", 
       y = "Marginal effect of SOLS change",
       title = "International constraint PCA score",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect. Darker shaded area",
                        " corresponds to estimated effect \nfrom minimum ",
                        " PCA score to ",
                        "2 standard deviations above mean."))

## Save plot
ggsave("appendix-figures/appx-sols-pca.pdf",
       height = 4,
       width = 8)


stargazer(pca_fit_1,
          dep.var.labels = c("Change in UNGA Ideal Point"),
          omit = c("ccode"),
          font.size = "tiny",
          covariate.labels = c(
            "SOLS Change",
            "Constraint score",
            "Other Leader Transition",
            "Democracy",
            "GDP per capita",
            "Population (log)",
            "Regime Transition",
            "Cold War End",
            "US Ally",
            "USSR Ally",
            "SOLS Change $\\times$ Constraint score",
            "Other Leader Transition $\\times$ Constraint score"
          ),
          label = "appx-pca-table",
          title = c(
            paste0("First principal component (constraint score)")),
          single.row = FALSE,
          no.space = TRUE,
          notes.align = "l",
          notes.label = "",
          df = FALSE,
          notes = c(
            "Two-tailed tests. Estimated standard errors in parentheses.", 
            "OLS estimates. Country dummies included in all models."),
          out = c("appendix-tables/appx-pca-table.tex",
                  "appendix-tables/appx-pca-table.html"))


#-----------------------------------------------------------------------
# Online Appendix: Accounting for interactions simultaneously
#-----------------------------------------------------------------------

all_fit <- lm(logvotech ~
                solsch * n_defenders +
                solsch * n_rivals +
                solsch * total_igos_impute +
                solsch * ln_total_budget_support +
                solsch * fdi_pct_gdp +
                solsch * trade_pct_gdp +
                solsch * state_capacity +
                othldrtrans * n_defenders +
                othldrtrans * n_rivals +
                othldrtrans * total_igos_impute +
                othldrtrans * ln_total_budget_support +
                othldrtrans * fdi_pct_gdp +
                othldrtrans * trade_pct_gdp +
                othldrtrans * state_capacity +
                demboth +
                gdp_per_capita +
                ln_population +
                regtrans + 
                CWend +
                USally +
                USSRally +
                factor(ccode) - 1, 
              data = dat[dat$interim == 0, ])

## Export to table
stargazer(all_fit,
          dep.var.labels = c("Change in UNGA Ideal Point"),
          omit = c("ccode", "othldrtrans"),
          font.size = "tiny",
          covariate.labels = c(
            "SOLS Change",
            "Defensive allies",
            "Rivals",
            "IGOs",
            "Aid",
            "FDI",
            "Trade",
            "Capacity",
            "Democracy",
            "GDP per capita",
            "Population (log)",
            "Regime Transition",
            "Cold War End",
            "US Ally",
            "USSR Ally",
            "SOLS Change $\\times$ Defensive allies",
            "SOLS Change $\\times$ Rivals",
            "SOLS Change $\\times$ IGOs",
            "SOLS Change $\\times$ Aid",
            "SOLS Change $\\times$ FDI",
            "SOLS Change $\\times$ Trade",
            "SOLS Change $\\times$ Capacity",
            "SOLS Change $\\times$ Democracy"
          )
          ,
          label = "appx-all-at-once-table",
          title = c(
            paste0("Including conditioning factors simultaneously")),
          
          single.row = FALSE,
          no.space = TRUE,
          notes.align = "l",
          notes.label = "",
          df = FALSE,
          notes = c(
            "Two-tailed tests. Estimated standard errors in parentheses.", 
            "OLS estimates. Country dummies included in all models."),
          out = c("appendix-tables/appx-all-at-once.tex", 
                  "appendix-tables/appx-all-at-once.html"))

## Get variable means
var_means <- model.frame(all_fit) %>%
  data.frame() %>%
  summarize_all(funs(mean(., na.rm = T))) %>%
  dplyr::select(n_defenders,
                n_rivals,
                total_igos_impute,
                ln_total_budget_support,
                fdi_pct_gdp,
                trade_pct_gdp,
                state_capacity)

## Get unique values from model
n_defenders_values <- unique(model.frame(all_fit)$n_defenders)
n_rivals_values <- unique(model.frame(all_fit)$n_rivals)
total_igos_impute_values <- unique(model.frame(all_fit)$total_igos_impute)
ln_total_budget_support_values <- unique(
  model.frame(all_fit)$ln_total_budget_support
)
fdi_pct_gdp_values <- unique(model.frame(all_fit)$fdi_pct_gdp)
trade_pct_gdp_values <- unique(model.frame(all_fit)$trade_pct_gdp)
state_capacity_values <- unique(model.frame(all_fit)$state_capacity)

## Put all in a list
all_u_vals <- list(
  n_defenders_values,
  n_rivals_values,
  total_igos_impute_values,
  ln_total_budget_support_values,
  fdi_pct_gdp_values,
  trade_pct_gdp_values,
  state_capacity_values
)

## Placeholder
calc_data <- list()

calc_data[[1]] <- data.frame(as.data.frame(var_means %>%
                                             dplyr::select(-n_defenders)), 
                             n_defenders = n_defenders_values,
                             model = "n_defenders")
calc_data[[2]] <- data.frame(as.data.frame(var_means %>%
                                             dplyr::select(-n_rivals)), 
                             n_rivals = n_rivals_values,
                             model = "n_rivals")
calc_data[[3]] <- data.frame(as.data.frame(var_means %>%
                                             dplyr::select(-total_igos_impute)), 
                             total_igos_impute = total_igos_impute_values,
                             model = "total_igos_impute")
calc_data[[4]] <- data.frame(
  as.data.frame(
    var_means %>%
      dplyr::select(-ln_total_budget_support)), 
  ln_total_budget_support = ln_total_budget_support_values,
  model = "ln_total_budget_support")
calc_data[[5]] <- data.frame(as.data.frame(var_means %>%
                                             dplyr::select(-fdi_pct_gdp)), 
                             fdi_pct_gdp = fdi_pct_gdp_values,
                             model = "fdi_pct_gdp")
calc_data[[6]] <- data.frame(as.data.frame(var_means %>%
                                             dplyr::select(-trade_pct_gdp)), 
                             trade_pct_gdp = trade_pct_gdp_values,
                             model = "trade_pct_gdp")
calc_data[[7]] <- data.frame(as.data.frame(var_means %>%
                                             dplyr::select(-state_capacity)), 
                             state_capacity = state_capacity_values,
                             model = "state_capacity")

calc_data <- bind_rows(calc_data)

## Putting them in order
for_mult <- calc_data %>%
  mutate(constant = 1) %>%
  dplyr::select(constant,
                n_defenders,
                n_rivals,
                total_igos_impute,
                ln_total_budget_support,
                fdi_pct_gdp,
                trade_pct_gdp,
                state_capacity)


library(MASS)
set.seed(32520)  ## Set seed to date

coef_sim <- cbind.data.frame(mvrnorm(n = 1000, 
                                     mu = coef(all_fit),
                                     Sigma = vcov(all_fit)) %>% 
                               as.data.frame() %>%
                               dplyr::select(contains("solsch")))

est_list <- list()

for(i in 1:nrow(coef_sim)){
  
  est_list[[i]] <- as.matrix(for_mult) %*% t(as.matrix(coef_sim[i, ]))
  
  
}


estimates <- bind_cols(est_list) %>% as.matrix()

calc_data$mean <- rowMeans(estimates)

calc_data$lower_95 <- sapply(1:nrow(estimates), 
                             function(x)
                               quantile(estimates[x, ] %>% as.numeric(), .025)[1])

calc_data$upper_95 <- sapply(1:nrow(estimates), 
                             function(x)
                               quantile(estimates[x, ] %>% as.numeric(), .975)[1])



## 1. Defensive allies
ggplot(calc_data %>% filter(model == "n_defenders"),
       aes(x = n_defenders,
           y = mean,
           ymin = lower_95,
           ymax = upper_95)) +
  geom_line(stat = "identity") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Defensive allies", 
       y = "Marginal effect of SOLS change",
       title = "Defensive alliances",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect."))

getwd()
ggsave(filename = "appendix-figures/appx-simultaneous-n-defenders.pdf",
       height = 4,
       width = 8)



## 2. Rivals
ggplot(calc_data %>% filter(model == "n_rivals"),
       aes(x = n_rivals,
           y = mean,
           ymin = lower_95,
           ymax = upper_95)) +
  geom_line(stat = "identity") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Rivalries", 
       y = "Marginal effect of SOLS change",
       title = "Rivalries",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect."))

getwd()
ggsave(filename = "appendix-figures/appx-simultaneous-n-rivals.pdf",
       height = 4,
       width = 8)


## 3. IGOs
ggplot(calc_data %>% filter(model == "total_igos_impute"),
       aes(x = total_igos_impute,
           y = mean,
           ymin = lower_95,
           ymax = upper_95)) +
  geom_line(stat = "identity") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Total IGO memberships", 
       y = "Marginal effect of SOLS change",
       title = "IGO memberships",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect."))

getwd()
ggsave(filename = "appendix-figures/appx-simultaneous-n-igos.pdf",
       height = 4,
       width = 8)


## 4. Aid
ggplot(calc_data %>% filter(model == "ln_total_budget_support"),
       aes(x = ln_total_budget_support,
           y = mean,
           ymin = lower_95,
           ymax = upper_95)) +
  geom_line(stat = "identity") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Total aid (budget support, log)", 
       y = "Marginal effect of SOLS change",
       title = "Foreign aid",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect."))

getwd()
ggsave(filename = "appendix-figures/appx-simultaneous-foreign-aid.pdf",
       height = 4,
       width = 8)



## 5. FDI
ggplot(calc_data %>% 
         filter(model == "fdi_pct_gdp") %>%
         filter(fdi_pct_gdp >= 0 & (fdi_pct_gdp <= 100)),
       aes(x = fdi_pct_gdp,
           y = mean,
           ymin = lower_95,
           ymax = upper_95)) +
  geom_line(stat = "identity") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "FDI as pct. GDP", 
       y = "Marginal effect of SOLS change",
       title = "FDI/GDP",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect."))

getwd()
ggsave(filename = "appendix-figures/appx-simultaneous-fdi.pdf",
       height = 4,
       width = 8)


## 6. Trade
ggplot(calc_data %>% filter(model == "trade_pct_gdp"),
       aes(x = trade_pct_gdp,
           y = mean,
           ymin = lower_95,
           ymax = upper_95)) +
  geom_line(stat = "identity") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Trade as pct. GDP", 
       y = "Marginal effect of SOLS change",
       title = "Trade/GDP",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect."))

getwd()
ggsave(filename = "appendix-figures/appx-simultaneous-trade.pdf",
       height = 4,
       width = 8)


## 7. Capacity
ggplot(calc_data %>% filter(model == "state_capacity"),
       aes(x = state_capacity,
           y = mean,
           ymin = lower_95,
           ymax = upper_95)) +
  geom_line(stat = "identity") +
  geom_hline(yintercept = 0, lty = "dashed", size = .25) +
  geom_ribbon(alpha = .3) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "State capacity", 
       y = "Marginal effect of SOLS change",
       title = "State capacity",
       caption = paste0("\nNote: Shaded region indicates 95% confidence ",
                        "interval around \nmarginal effect."))

getwd()
ggsave(filename = "appendix-figures/appx-simultaneous-capacity.pdf",
       height = 4,
       width = 8)


#-----------------------------------------------------------------------
# Online Appendix: 3-way interaction between SOLS, conditions, population
#-----------------------------------------------------------------------

## For figure in Section J of the Appendix



# Model 1: Number of defensive allies
sfit_3waypop_1 <- lm(logvotech ~
                       solsch * n_defenders * ln_population +
                       othldrtrans * n_defenders * ln_population + 
                       demboth + 
                       gdp_per_capita + 
                       regtrans + 
                       CWend +
                       USally +
                       USSRally +
                       factor(ccode) - 1, 
                     data = dat[dat$interim == 0, ])

## Get upper bound of number of defensive allies to consider
n_def_max <- round(mean(sfit_3waypop_1$model$n_defenders) +
                     2 * sd(sfit_3waypop_1$model$n_defenders), 0)

defenders_range <- c(0, n_def_max)

# Model 2: Number of rivals

sfit_3waypop_2 <- lm(logvotech ~
                       solsch * n_rivals * ln_population +
                       othldrtrans * n_rivals * ln_population + 
                       demboth + 
                       gdp_per_capita +
                       regtrans + 
                       CWend +
                       USally +
                       USSRally +
                       factor(ccode) - 1, 
                     data = dat[dat$interim == 0, ])

n_rivals_max <- round(mean(sfit_3waypop_2$model$n_rivals) +
                        2 * sd(sfit_3waypop_2$model$n_rivals), 0)

rivals_range <- c(0, n_rivals_max)

## State capacity -------------------------------------------------------------

capacity_fit_3waypop_1 <- lm(logvotech ~
                               solsch * state_capacity * ln_population +
                               othldrtrans * state_capacity * ln_population + 
                               demboth + 
                               gdp_per_capita +
                               regtrans + 
                               CWend +
                               USally +
                               USSRally +
                               factor(ccode) - 1, 
                             data = dat[dat$interim == 0, ])


capacity_max <- round(mean(capacity_fit_3waypop_1$model$state_capacity) +
                        2*sd(capacity_fit_3waypop_1$model$state_capacity), 0)

capacity_min <- round(mean(capacity_fit_3waypop_1$model$state_capacity) -
                        2*sd(capacity_fit_3waypop_1$model$state_capacity), 0)

capacity_range <- c(capacity_min, capacity_max)

## IGO model ----------------------------------------------------------------
igo_fit_3waypop_1 <- lm(logvotech ~
                          solsch * total_igos_impute * ln_population +
                          othldrtrans * total_igos_impute * ln_population + 
                          demboth + 
                          gdp_per_capita +
                          regtrans + 
                          CWend +
                          USally +
                          USSRally +
                          factor(ccode) - 1, 
                        data = dat[dat$interim == 0, ])

## Get upper bound for + 2 SDs
n_igos_max <- round(mean(igo_fit_3waypop_1$model$total_igos_impute) +
                      2 * sd(igo_fit_3waypop_1$model$total_igos_impute), 0)


igos_range <- c(0, n_igos_max)

## Political economy models -------------------------------------------------
aid_fit_3waypop_1 <- lm(logvotech ~
                          solsch * ln_total_budget_support * ln_population +
                          othldrtrans * ln_total_budget_support * ln_population + 
                          demboth + 
                          gdp_per_capita + 
                          regtrans + 
                          CWend +
                          USally +
                          USSRally +
                          factor(ccode) - 1, 
                        data = dat[dat$interim == 0, ])


## Get upper bound for + 2 SDs
aid_max <- round(mean(aid_fit_3waypop_1$model$ln_total_budget_support) +
                   2 * sd(aid_fit_3waypop_1$model$ln_total_budget_support), 0)

aid_min <- round(mean(aid_fit_3waypop_1$model$ln_total_budget_support) -
                   2 * sd(aid_fit_3waypop_1$model$ln_total_budget_support), 0)

## Minus 2 SDs is -11, replacing as 0
aid_min <- 0

aid_range <- c(aid_min, aid_max)



## Foreign direct investment ------------------------------------------------- 
fdi_fit_3waypop_1 <- lm(logvotech ~
                          solsch * fdi_pct_gdp * ln_population +
                          othldrtrans * fdi_pct_gdp * ln_population + 
                          demboth + 
                          gdp_per_capita + 
                          regtrans + 
                          CWend +
                          USally +
                          USSRally +
                          factor(ccode) - 1, 
                        data = dat[dat$interim == 0, ])


fdi_max <- round(mean(fdi_fit_3waypop_1$model$fdi_pct_gdp) +
                   2*sd(fdi_fit_3waypop_1$model$fdi_pct_gdp), 0)

fdi_min <- round(mean(fdi_fit_3waypop_1$model$fdi_pct_gdp) -
                   2*sd(fdi_fit_3waypop_1$model$fdi_pct_gdp), 0)

fdi_range <- c(0, fdi_max)


## Trade dependence -----------------------------------------------------------

trade_fit_3waypop_1 <- lm(logvotech ~
                            solsch * trade_pct_gdp * ln_population +
                            othldrtrans * trade_pct_gdp * ln_population + 
                            demboth + 
                            gdp_per_capita + 
                            regtrans + 
                            CWend +
                            USally +
                            USSRally +
                            factor(ccode) - 1, 
                          data = dat[dat$interim == 0, ])


trade_max <- round(mean(trade_fit_3waypop_1$model$trade_pct_gdp) +
                     2*sd(trade_fit_3waypop_1$model$trade_pct_gdp), 0)

trade_min <- round(mean(trade_fit_3waypop_1$model$trade_pct_gdp) -
                     2*sd(trade_fit_3waypop_1$model$trade_pct_gdp), 0)

trade_range <- c(trade_min, trade_max)




#-----------------------------------------------------------------------
# Figures
#-----------------------------------------------------------------------

ranges <- list(defenders_range,
               rivals_range,
               capacity_range,
               igos_range,
               trade_range,
               aid_range,
               fdi_range)

fits <- list(sfit_3waypop_1,
             sfit_3waypop_2,
             capacity_fit_3waypop_1,
             igo_fit_3waypop_1,
             trade_fit_3waypop_1,
             aid_fit_3waypop_1,
             fdi_fit_3waypop_1)


conditioning_values <- list(
  seq(defenders_range[1], defenders_range[2], 
      length.out = max(c(20, defenders_range[2]))),
  seq(rivals_range[1], rivals_range[2], 
      length.out = max(c(20, rivals_range[2]))),
  seq(capacity_range[1], capacity_range[2], 
      length.out = max(c(20, capacity_range[2]))),
  seq(igos_range[1], igos_range[2], 
      length.out = max(c(20, igos_range[2]))),
  seq(trade_range[1], trade_range[2], 
      length.out = max(c(20, trade_range[2]))),
  seq(aid_range[1], aid_range[2], 
      length.out = max(c(20, aid_range[2]))),
  seq(fdi_range[1], fdi_range[2], 
      length.out = max(c(20, fdi_range[2])))
)


titles <- c("Number of defensive allies",
            "Number of rivals",
            "State capacity",
            "Number of IGO memberships",
            "Trade",
            "Aid dependence",
            "FDI")

pop_upper <- round(mean(dat$ln_population, na.rm = T) +
                     2*sd(dat$ln_population, na.rm = T), 2)

pop_range <- c(min(dat$ln_population, na.rm = T), pop_upper)

pop_range_values <- seq(pop_range[1], pop_range[2], length.out = 30)


## Simulate betas with mean = coefs from model and model var-cov matrix
library(MASS)

separate_figures <- list()

for(i in 1:length(fits)){
  
  sim_beta <- mvrnorm(n = 1000, 
                      mu = coef(fits[[i]]),
                      Sigma = vcov(fits[[i]]))
  
  ## Keep necessary columns
  sim_beta <- sim_beta %>%
    as.data.frame() %>%
    dplyr::select(contains("solsch"))
  
  ## Place holder for each data set
  calc_data <- list()
  
  for(j in 1:nrow(sim_beta)){
    
    betas <- sim_beta[j, ] %>% unlist() %>% c()
    
    ## Set up data set to look at all combos of condition / democracy
    calc_data[[j]] <- cbind.data.frame(
      condition = rep(conditioning_values[[i]], length(pop_range_values)),
      Population = c(rep(pop_range_values,  each = length(conditioning_values[[i]]))),
      b_1 = betas[1],
      b_2 = betas[2],
      b_3 = betas[3],
      b_4 = betas[4],
      row.names = NULL
    ) %>%
      mutate(
        beta_sols = b_1 + 
          (b_2 * condition) + 
          (b_3 * Population) + 
          (b_4 * condition * Population)
      ) %>%
      dplyr::select(condition, Population, beta_sols) %>%
      mutate(iteration = j)
    
  }
  
  ## Bind all together
  calc_data <- do.call(rbind, calc_data) %>%
    as.data.frame()
  
  ## Calculate mean and lower/upper quantile
  calc_data <- calc_data %>%
    group_by(condition, Population) %>%
    summarise(mean_beta = mean(beta_sols),
              beta_lower = quantile(beta_sols, 0.025),
              beta_upper = quantile(beta_sols, 0.975)) %>%
    ungroup() %>%
    mutate(significant = ifelse(beta_lower * beta_upper > 0, TRUE, FALSE)) %>%
    mutate(positive = ifelse(mean_beta > 0, TRUE, FALSE)) %>%
    mutate(Result = NA) %>%
    mutate(Result = replace(Result,
                            positive & significant,
                            "+, significant")) %>%
    mutate(Result = replace(Result,
                            positive & (!significant),
                            "+")) %>%
    mutate(Result = replace(Result,
                            (!positive) & (!significant),
                            "-")) %>%
    mutate(Result = replace(Result,
                            (!positive) & significant,
                            "-, significant")) %>%
    mutate(Result = factor(Result,
                           levels = rev(c("-, significant",
                                          "-",
                                          "+",
                                          "+, significant")))) 
  
  
  
  p <- ggplot(calc_data,
              aes(x = condition, y = Population,
                  fill = Result)) +
    scale_fill_manual(values = c("#383838",
                                 "#808080",
                                 "#d3d3d3",
                                 "#f5f5f5")) +
    geom_tile() +
    theme_plot_plain() +
    labs(
      x = titles[i],
      y = "Population (log)"
    ) 
  
  ## Uncomment to save each panel separately
  # getwd()
  # the_file <- paste0("figures/",
  #                    titles[i] %>% tolower() %>% gsub(" ", "-", .),
  #                    "-by-population",
  #                    ".jpeg")
  # 
  # 
  # ## Save
  # ggsave(p,
  #        file = the_file)
  
  separate_figures[[i]] <- p
  
}


## Make full plot
library(gridExtra)

ggsave(
  plot = grid.arrange(
    separate_figures[[1]],
    separate_figures[[2]],
    separate_figures[[3]],
    separate_figures[[4]],
    separate_figures[[5]],
    separate_figures[[6]],
    separate_figures[[7]], 
    ncol = 2
  ),
  filename = "appendix-figures/pop-3way-figure.png",
  width = 12, height = 16
)


